Phase behavior and structure of colloidal bowl-shaped particles: simulations 
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We study the phase behavior of bowl-shaped particles using computer simulations. These particles 
were found experimentally to form a meta-stable worm-like fluid phase in which the bowl-shaped 
particles have a strong tendency to stack on top of each other [M. Marechal et al, Nano Letters 
10, 1907 (2010)]. In this work, we show that the transition from the low-density fluid to the 
worm-like phase has an interesting effect on the equation of state. The simulation results also 
show that the worm-like fluid phase transforms spontaneously into a columnar phase for bowls that 
are sufficiently deep. Furthermore, we describe the phase behavior as obtained from free energy 
calculations employing Monte Carlo simulations. The columnar phase is stable for bowl shapes 
ranging from inflnitely thin bowls to surprisingly shallow bowls. Aside from a large region of 
stability for the columnar phase, the phase diagram features four novel crystal phases and a region 
where the stable fluid contains worm-like stacks. 
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I. INTRODUCTION 

The concept of a mesogenic particle in the form of a 
bowl is relatively old in the molecular liquid crystal com- 
munity. Such molecules are expected to form a columnar 
phase, which can be ferroelectric, i.e., a phase with a 
net electric dipolc moment, when the particles possess a 
permanent dipole moment. Ferroelectric phases have po- 
tential applications for optical and electronic devices. In 
fact, crystalline (as opposed to liquid crystalline) ferro- 
electrics are already applied in sensors, electromechanical 
devices and non- volatile memory jlj. A columnar ferro- 
electric phase may have the advantage over a crystal, that 
grain boundaries and other defects anneal out faster due 
to the partially fluid nature of the columnar phase. In re- 
ality, columnar phases of conventional disc-like particles 
often exhibit many defects, as flat thin discs can diffuse 
out of a column and columns can split up. The presence 
of these defects limits their potential use for industrial 
applications [2]. Less defects are expected in a columnar 
phase of bowl-shaped mesogens, where particles are sup- 
posed to be more confined in the lateral directions. A 
whole variety of bowl-like molecules have already been 
synthesized and investigated experimentally [MS]. In 
addition, buckybowlic molecules, i.e. fragments of Cgo 
whose dangling bonds have been saturated with hydro- 
gen atoms, have been shown to crystallize in a columnar 
fashion [7HTT] . However, the number of theoretical stud- 
ies is very limited as it is difficult to model the compli- 
cated particle shape in theory and simulations. In a re- 
cent simulation study, the attractive-repulsive Gay-Berne 
potential generalized to bowl-shaped particles has been 
used to investigate the stacking of bowl-like mesogens as 
a function of temperature |2j . The authors reported a ne- 
matic phase and a columnar phase. This columnar phase 
did not exhibit overall ferroelectric order, although polar 
regions were found. In another very recent simulation 
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study [12] of hard contact lenses (infinitely thin, shallow 
bowls), a new type of fluid phase was found in which the 
particles cluster on a spherical surface for bowls which 
are not too shallow. No columnar phase was found since 
the focus was on rather shallow bowls at a relatively low 
densities. 

Recently, a procedure has been developed to synthesize 
bowl-shaped colloidal particles [13]. This method starts 
with the preparation of highly uniform oil-in-water emul- 
sion droplets. Subsequently, the droplets were used as 
templates around which a solid shell with tunable thick- 
ness is grown. In the next step of the synthesis, the oil in 
the droplets is dissolved and finally, during drying, the 
shells collapse into hemispherical double- walled bowls. In 
addition to these larger, more easily imaged colloids, a 
whole variety of bowl-shaped nanoparticles and smaller 
colloids have been synthesized and characterized [TIHT^ , 
and possible applications of these systems have been put 
forward. We also note that recently hemispherical parti- 
cles were synthesized at an air-solution interface [20] and 
on a substrate [H]. These hemispherical particles are 
intended to be used as microlense arrays, but they can 
also serve as a new type of shape-anisotropic colloidal 
particle. 

In our simulations, we model the particles as the solid 
of revolution of a crescent (see Fig. [l^). The diameter 
a of the particle and the thickness D are defined as in- 
dicated in Fig. [T^. We define the shape parameter of 
the bowls by a reduced thickness D/a, such that the 
model reduces to infinitely thin hemispherical surfaces 
for D /a = and to solid hemispheres for D/a = 0.5. 
The advantages of this simple model is that it interpo- 
lates continuously between an infinitely thin bowl and 
a hemispherical solid particle (the two colloidal model 
systems, which we discussed above), and that we can de- 
rive an algorithm that tests for overlaps between pairs 
of bowls, which is a prerequisite for Monte Carlo simula- 
tions of hard-core systems. 

In a recent combined experimental and simulation 
study (for which we performed the simulations), the 




FIG. 1. (a) The theoretical model of the colloidal bowl is 
the solid of revolution of a crescent around the axis indicated 
by the dashed line. The thickness of the double- walled bowl 
is denoted by D and the diameter of the bowl by a. (b) The 
bowls are defined using two spheres of radii -Ri and R2 , that 
are a distance of L apart. The direction vector, Ui and the 
reference point of the particle, r^, (the dot in the center of the 
smaller sphere) are indicated. 



phase behavior of repulsive bowl-shaped colloids was in- 
vestigated [H] . The colloids were shown to form a worm- 
hke fluid phase, in which the particles form long curved 
stacks running in random directions. By comparing the 
distribution of stack lengths, the simulation model was 
shown to describe the colloidal particles well. No evi- 
dence of columnar ordering was found in the experiments 
and in simulations of bowls with corresponding deepness, 
which was explained by the glassy behavior of the par- 
ticles preventing rearrangements. The phase behavior of 
the model particles is expected to also describe other re- 
pulsive bowl shaped particles well, provided that the di- 
mensions of the simulation particle are chosen such that 
the diameter of a stack and the inter-particle distance in 
the stack are the same as for the particles to be modeled. 
In this work, we expand the simulation results on the 
hard bowl-shaped particles. First, we elaborate on the 
model for the collapsed shells; the overlap algorithm is 
described in the appendix. Also, the (free energy) meth- 
ods are explained in more detail than in Ref. 22 . In the 
results section, we study the properties of the isotropic 
phase. We investigate the nature and the location of the 
transition between the homogeneous fluid phase and the 
fluid phase that contains the worm-like stacks. Further- 
more, we show the packing diagram and the phase di- 
agram with a tentative homogeneous-to-worm-like fluid 
transition line. In the last section we summarize and 
discuss the results. 



II. METHODS 

A. Model 

We describe the model that we use to represent the 
bowls in more detail. Consider a sphere with a radius Ri 
at the origin and a second sphere with radius R2 > Ri 
at position —Lui, where u^ is the unit vector denoting 
the orientation of the bowl and L > 0. The bowl is 
represented by that part of the sphere with radius Ri that 
has no overlap with the larger sphere, see Fig. ^p. We 



have chosen the values for L and R2 such that the bowls 
are hemispherical (see appendix for explicit expressions 
for L and i?2)- We define the thickness of the bowls 
hy D — L — (i?2 — ^1), such that the model reduces to 
the surface of a hemisphere for D = and to a solid 
hemisphere for D = Ri. The volume of the particle is 
jD{a'^ - Da + ^D^), where a = 2Ri is our unit of 
length. The algorithm to determine overlap between our 
bowls is described in the appendix. 



B. Fluid phase 

We employ standard NPT MC simulations to obtain 
the equation of state (EOS) for the fluid phase. In ad- 
dition, we obtain the compressibility by measuring the 
fluctuations in the volume: 



{V) 



P 



dp 
dP 



(1) 



where p = N/V is the number density and the deriva- 
tive of the pressure is taken at constant temperature is 
denoted by the subscript T. We determine the free en- 
ergy at density pi by integrating the EOS from reference 
density po to pi: 
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where the chemical potential /i(/9o) is determined using 
the Widom particle insertion method [23], and P{po) is 
determined by a local fit to the EOS. 

To investigate the structure of the fluid phase, we mea- 
sure the positional correlation function [24| . 



9c{z) 



1 



NpA, 



N WoolCO 



Uj-z)), 



(3) 



where the sum over j runs over Nco\{i) particles in a col- 
umn of radius ct/2 with orientation Ui centered around 
particle i, and where the area of the column is denoted by 
^coi = 7rcr^/4. At sufficiently high pressure the particles 
stack on top of each other to form disordered worm-like 
piles which resemble the stacks observed in the exper- 
iments [22 . As the stacks have a strong tendency to 
buckle, we cannot use gdz) to determine the length of 
the stacks. We therefore determine the stack size distri- 
bution using a cluster criterion. Particle i and j belong 
to the same cluster if 



r„- + {CD/2 + cr/4)(uj - Uj)| < cr/2 and 

U,; • Uj > 0, (4) 



and where the first condition has to be satisfied for ^ = 
— 1, or 1 and r^ = Tj — Vi, with r^ denoting the center 
of the sphere with radius i?i of particle i, see Fig. [I]d. 
If both conditions are satisfied, particle j is just above 
(C = 1) or below {( = —1) particle i in the stack, or, when 



the stack is curved, particle j can be next to particle i 
(C = 0). We now define the cluster distribution as the 
fraction of particles that belongs to a cluster of size n: 
Vstackin) = nNn/N, where Nn is the number of clusters 
of size n. We checked that the cluster size distribution 
does not depend sensitively to the choice of parameters 
inEq. ^. 



C. Columnar phases 

We also perform NPT Monte Carlo simulations of the 
columnar phase using a rectangular simulation box with 
varying box lengths in order to relax the inter-particle 
distance in the z direction, along the columns, indepen- 
dently from the lattice constant in the horizontal direc- 
tion. The difference between the free energy of the colum- 
nar phase at a certain density and the free energy of 
the fluid phase at a lower density is determined using a 
thermodynamic integration technique [5S]- We apply a 
potential which couples a particle to its column: 
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where Xi and yi are the x and y components respectively 
of Ti, Na is the number of columns in the a direction 
and La is the size of the box in the a direction. In our 
simulations, we calculate Eq. (Is]) while fixing the center 
of mass. To do so efficiently, we first calculate all four 
combinations 
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for trigl = cos, sin and trig2 = cos, sin. The change 
in these four expressions upon displacement of a single 
particle while keeping the center of mass fixed can be 
expressed in terms of single particle properties and the 
previous values of the expressions by using some basic 
trigonometry. In this way, <i>hcx(r^. A), which is Eq. [6] 
for trigl — cos and trig2 = sin, can be calculated with- 
out performing the full summation over all particles in 
Eq. (|6| every time we displace a particle. Unfortunately, 
this calculation requires the evaluation of many more 
trigonometric functions than the simple expression ([5]), 
but the extra computation time is negligible compared 
to the overlap check. 

In addition to this positional potential, we also con- 
strain the direction of the particle, using the potential 
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from the columnar phase at a certain density p2- Subse- 
quently, we slowly turn on the two potentials, i.e. we in- 
crease A from to Amax- Next, we integrate the equation 
of state to go from p2 to pi , while keeping A = Amax fixed. 
During this step the columnar phase will only be stable 
below the coexistence density, if Amax is sufficiently high. 
We find that A^ax = ^OkgT suffices to guarantee stability 
of the columnar phase. Finally, fixing the density pi , we 
gradually turn off the potentials, while integrating over 
A from Amax to 0. During this last step, the columnar 
phase melts continuously, provided that the density pi is 
low enough and that A is high enough to prevent melting 
during the density integration step. The resulting free 
energy difference between the columnar phase and fluid 
phase is given by 

Fcol{p2) --Ffluid(/5l) = 
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The positional potential (l5| is designed to stabilize a 
hexagonal array of columns, but, strictly speaking, it 
does not have the hexagonal symmetry of the columnar 
phase, since it is not invariant under a 60 degrees rotation 
of the whole system around a lattice position. However, 
we find that replacing Eq. ([5]) by a positional potential 
that does have this symmetry, does not have a significant 
effect on the free energy difference. 

A second type of columnar phase can be constructed 
by flipping half of the bowls. In this way we obtain al- 
ternating vertical sheets (i.e. rows of columns) of bowls 
that point upwards and sheets of bowls that point down- 
wards, we will refer to this phase as the inverted columnar 
phase. We calculate the free energy of this phase using 
the method described above, with the modification that 
the angular potential now reads, 



*angK,A) = A' 
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This potential could also have been used for the non- 
inverted columnar phase, and we have found that the re- 
sult of the free energy integration for the columnar phase 
is the same whether we use Eq. ^ or Eq. ([t]). 



D. Crystals 



1. Packing 



where we used A' = 0.1 A and where Ui^z is the z com- 
ponent of Uj. The thermodynamic integration path from 
the columnar phase to the fluid is as follows: We start 



As the crystal phases of the bowls are not known a pri- 
ori, we developed a novel pressure annealing method to 
obtain the possible crystal phases |26j, which we named 



after the thermal anneaUng technique commonly used 
to find energy minima. Fully variable box shape NPT 
simulations were performed on system of only 2-6 par- 
ticles. By construction, the final configuration of such 
a simulation is a crystal, where the unit cell is the 
simulation box. One cycle of such a simulation con- 
sists of the following steps: We start at a pressure of 
lOkBT/a^. Subsequently, we run a series of simulations, 
where the pressure increases by a factor of ten each run: 
Pa^/keT = 10, 100, . . . , 10^ At the highest pressure 
(lO^ksT/a^) we measure the density and angular order 
parameters, Si = ||(u,i)|| and 5*2 = A2, where A2 is the 
highest eigenvalue of the matrix whose components are 
Qap = ^{uiaUip) — \5ap, whcrc a, /3 = x,y,z. We store 
the density if it is the highest density found so far for 
these values of 5i and S'2. We ran 1000 of such cycles for 
each aspect ratio, which is enough to visit each crystal 
phase multiple times. After completing the simulations, 
we tried to determine the lattice parameters of the re- 
sulting crystal by hand. Although this last step is not 
necessary, it is convenient to have analytical expressions 
for the lattice vectors and the density. The pressure an- 
nealing runs were performed for D/a = 0.1, 0.15, . . . , 0.5. 
For many of the crystals, we were not able to find ana- 
lytical expressions for the lattice parameters. For these 
crystals, we obtain the densities of the close packed crys- 
tals for intermittent values of D by averaging the density 
in single simulation runs at a pressure of lO^ksT/a^. 
The initial configurations for the value of D of interest 
were obtained from the final configurations of the pres- 
sure annealing simulations for another value of D by one 
of the following two methods, depending on whether we 
needed to decrease or increase D: When decreasing L no 
overlaps are created so the final configuration of the sim- 
ulation for the previous value of L can be used as initial 
configuration. On the other hand, increasing L results in 
an overlap, which is removed by scaling the system uni- 
formly. Subsequently, the pressure is stepwise increased 
from lOOOksT/a^ to lO^ksT/a^, by multiplying by 10 
each step. 



tential by the interaction strength 7. We integrate over 
dF/dj from a system with essentially hard-core interac- 
tion (high 7 = 7max), to an ideal Einstein crystal (7 = 0). 
Some minor alterations to the scheme of Ref. [55] were in- 
troduced, which were necessary, because of the different 
shape of the particle. For the coupling of the orientation 
of bowl i, i.e., Uj, to an aligning field, we have to take into 
account that the bowls have no up down symmetry, while 
the dumbbells are symmetric under u^ — ?► — u^. The po- 
tential energy function that achieves the usual harmonic 
coupling of the particles to their lattice positions, as well 
as the new angular coupling, reads: 
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where r^ and u^ denote, respectively, the center-of-mass 
position and orientation of bowl i and ro,i the lattice site 
of particle i, 9io is the angle between u^ and the ideal 
tilt vector of particle i, and f3 — l/ksT. The Hclmholtz 
free energy |28| of the noninteracting Einstein crystal is 
modified accordingly, but the only modification is the 
integral over the angular coordinates: 
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Although the shape of the bowls is more complex than 
that of the dumbbell, we can still use a rather simple 
form for the pairwise soft potential interaction: 

/3C/,oft(r^,u^;7) = EMr,-r„u„u„7) (12) 



i<j 



with 



/3(^(rj -r,,Ui,Uj,7) = 

7(1 - A(r^ /(Tinax)^) if i and j overlap 
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2. Free energies 

We calculate the free energy of the various crystal 
phases by thermodynamic integration using the Einstein 
crystal as a reference state ^7^. The Einstein integra- 
tion scheme that we employ here is similar to the one 
that was used to calculate the free energies of crystals 
of dumbbells in Ref. ^8, . We briefiy sketch the integra- 
tion scheme here and discuss the modifications that we 
applied. We couple both the positions and the direction 
of the particles with a coupling strength A, such that for 
A — )■ 00, the particles are in a perfect crystalline configu- 
ration. First, we integrate dF/dX over A from zero to a 
large but finite value for A. Subsequently we replace the 
hard-core particle-particle interaction potential by a soft 
interaction, where we can tune the softness of the po- 



where r^ = |rj - r^ -(- ^-^ (u^ - u 
tween the "centers" of bowl 
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i.e. the distance be- 
i and bowl j, CTmax is the 
for which the particles overlap: cTmax ^ 
a'"^ + {a — DY , A is an adjustable parameter that is kept 
fixed during the simulation at a value A = 0.5, and 7 
is the integration parameter. It was shown in Ref. |29j 
that in order to minimize the error and maximize the ef- 
ficiency of the free energy calculation, the potential must 
decrease as a function of r and must exhibit a disconti- 
nuity at r such that both the amount of overlap and the 
number of overlaps decrease upon increasing 7. Here, 
we have chosen this particular form of the potential be- 
cause it can be evaluated very efficiently in a simulation, 
although it does not describe the amount of overlap be- 
tween bowls i and j very accurately. We checked that 
adding a term that tries to describe the angular behavior 




FIG. 2. The final configuration obtained from simulations 
at Pa^ /ksT = 50 and D = O.Sa The colors denote different 
stacks. 



of the amount of overlap does not significantly change our 
results of the free energy calculations. Also, we checked 
that by employing the usual Einstein integration method 
(i.e. only hard-core interactions) at a relatively low den- 
sity we obtained the same result as by using the method 
of Fortini et al. |29j . Finally, we set the maximum inter- 
action strength 7i„ax to 200. 

We perform variable box shape NPT simulations [30] 
to obtain the equation of state for varying D. In these 
simulations not only the edge length changes, but also 
the angles between the edges are allowed to change. We 
employ the averaged configurations in the Einstein crys- 
tal thermodynamic integration. We calculate the free 
energy as a function of density by integrating the EOS 
from a reference density to the density of interest: 



F(pl) = ^(po) 



y{^) ,., 



III. RESULTS 



A. Stacks 



tonically with packing fraction (j) for D = O.lcr, see 
Fig. |3[ where the packing fraction is defined as = 
^((7^- Do- + lD'^)N/V. This behavior persist for all 
D < 0.2a, but for D > 0.25cr the pressure is always con- 
vex. We investigate the origin of these peculiarities using 
gdz), the positional correlation function along the direc- 
tor of a particle, which includes only the particles in a 
column around a particle, as defined in Eq. pi). As can 
be seen from gdz) in Fig. l4J the structure of the fluid 
changes dramatically as the pressure is increased. At 
P* = f3Pa^ = 1, the correlation function is typical for a 
low density isotropic fluid of hemispherical particles; no 
effect of the dent of the particles is found at low densities. 
The only peculiar feature of gdz) for P* = 1 is that it 
is not symmetric around zero, but this is caused by our 
choice of reference point on the particle (see Fig. IIId), 
which is located below the particle if the particle points 
upwards. In contrast, at P* — 10 gd^) already shows 
strong structural correlations. Most noteworthy is the 
peak at z = D, that shows that the fluid is forming short 
stacks of aligned particles. Also, note that the value of 
gd^) is nonzero around z = 0. This is caused by pairs 
of bowls that align anti-parallel and form a sphere-like 
object, as depicted in Fig. [4] Finally, at P* = 50 and 
higher, long worm-like stacks are fully formed and gdz) 
shows multiple peaks at z = Dn for both positive and 
negative integer values of n. Furthermore, at these pres- 
sures, there are no sphere-like pairs, as can be observed 
from the value of 5c (0). The formation of stacks explains 
the peculiar behavior of the pressure: At low densities, 
the bowls rotate freely, which means that the pressure 
will be dominated by the rotationally averaged excluded 
volume. The excluded volume of two particles that are 
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We perform standard Monte Carlo simulations in the 
isobaric- isothermal ensemble (NPT). Fig. [2] shows a typi- 
cal configuration of bowl-shaped particles with D — 0.3 a 
at Pa^ /ksT = 50, displaying stacking behavior typical 
for the worm-like phase. The equation of state (EOS) 
of the fluid is somewhat peculiar: the pressure as a 
function of density is not always convex for all densi- 
ties, although the compressibility does decrease mono- 



FIG. 3. The equation of state for bowl-shape particles with 
D — O.lcr, reduced pressure P* — pPa^ (left axis), and the 
reduced compressibility - Jj^ on a log scale (right axis) as a 
function of packing fraction (j). The points are data obtained 
from NPT simulations. The solid line is a fit to the pressure; 
the dashed line is the corresponding reduced compressibility, 

W 8fit(p) N-l 

p\ dp J 



not aligned is nonzero, even for D = 0, and gives rise to 
the convex pressure which is typical for repulsive parti- 
cles. As the density increases and the bowls start to form 
stacks, the available volume increases, and the pressure 
increases less than expected, which can even cause the 
EOS to be concave. At even higher densities the worm- 
like stacks are fully formed, and the pressure is again a 
convex function of density for D > 0, dominated by the 
excluded volume of locally aligned bowls. The excluded 
volume of completely aligned infinitely thin bowls is zero, 
and, therefore, the pressure increases almost linearly with 
density for D — when the stacks are fully formed. 

To quantify the length of the stacks we calculated the 
stack distribution as shown in Fig. [5) As can be seen 
from the figure, the length of the stacks is strongly de- 
pendent on D/a. However, we have found that above 
a certain threshold pressure the distribution of stacks is 
nearly independent of pressure. 

We investigated whether the worm-like stacks could 
spontaneously reorient to form a columnar phase. We 
increased the pressure in small steps of 1 kBT/a^ from 
well below the fluid-columnar transition to very high 
pressures, where the system was essentially jammed. At 
each pressure, we ran the simulation for 4 • 10^ Monte 
Carlo cycles, where a cycle consists of N particle and 
volume moves. These simulations show that the bowls 
with a thickness D > 0.25a always remained arrested in 
the worm-like phase, which is similar to the experimental 
observations |22|. However, for D/a — 0.1 and 0.2, we 
find that the system eventually transforms into a colum- 
nar phase in the simulations (see Fig. [6]) . This might be 




FIG. 4. The pair correlation function, gc{z), of a fluid of 
bowl-shaped particles with D = 0.2a as a function of the di- 
mensionless inter-particle distance z/a along the axis of a ref- 
erence bowl for various reduced pressures P* = 13 Pa^. Only 
particles within a cylinder of diameter a around the bowl are 
considered, as indicated by the subscript 'c'. We show typ- 
ical two-particle configurations that contribute to gc{z) for 
z/a = —0.5,-0.2,0.2,0.4 and 1, where the filled bowls de- 
note the reference particle, and the open bowls with thick 
outlines denote the other particle. 



o 




n 



FIG. 5. The probability, Pstmckin), to find a particle in a stack 
of size n for D/a = 0.2, 0.3 and 0.4 and Pa^/ksT = 50 



explained by the fact that the isotropic-to-columnar tran- 
sition occurs at lower packing fractions for deeper bowls 
(smaller D), which facilitates the rearrangements of the 
particles into stacks and the alignment of the stacks into 
the columnar phase. 



B. Packing 

We found six candidate crystal structures, denoted 
X, IX, IX', B, IB and fcc^, using the pressure annealing 
method. Snapshots of a few unit cells of these crystal 
phases are shown in Fig. [7] We will describe these crys- 




FIG. 6. The final configuration of a simulation of bowls with 
L = O.ID at Pa^/ksT = 38. The gray values denote diflterent 
columns. 





fee 



€» 



IX' 



FIG. 7. The various crystal phases that were considered as 
possible stable structures. Five of these were found using 
the pressure annealing method: X, IX, B, IB and IX'. X, 
IX, B and IB are densely packed structures for D < 0.5 a and 
fcc^ and IX' are densely packed crystal structures for (nearly) 
hemispherical bowls [D ~ 0.5 cr). 



tal structures using the order parameters 5*1, that mea- 
sures alignment of the particles, and the nematic order 
parameter (S'2), that is nonzero for both parallel and anti- 
parallel configurations. Crystal structure X has S*! ~ 1 
and 5'2 — 1, and the particles are stacked head to toe in 
columns. The lattice vectors are 



ai — ax 



3.2 — Dz 



-X + -\/cr2 _ i:»2 _^ 2a^fa^ - D^ 



ag = -a;+ -y a- - U- + 'Zax/ a- - U- y+ —z, 
and the density is 

'Da 

pa- 



a^-D^ + 2(7 Vcr2 _ £,2 



(15) 



(16) 



The order parameters of the second crystal structure, 
are Si ~ and S'2 — 1, which is caused by the fact 
that half of the particles point upwards, and the other 
half downwards. Further investigation shows that there 
are two phases with 5*2 — 1 and Si ~ 0: one at low D 
(IX) and one at _D ~ a/2 (IX'). The structure within 
the columns of the first (IX) of these two structures is 
the same as for the X structure, but one half of these 
columns are upside down, like in the inverted colum- 
nar phase (in fact, the IX crystal melts into the inverted 
columnar phase). The lattice vectors of crystal structure 
IX are 



ai = ax a2 = Dz 

cr . 1 



33 = -X + -Vi(r^ - ^D^ y, 



and the density is 



P^' = 



Da 



v/3ct2 - 4:D^ 



(17) 



(18) 



The columns in the IX crystal are arranged in such a way 
that the rims of the bowls can interdigitate. The IX' crys- 
tal can be obtained from the IX phase at D = a/2 by 
shifting every other layer by some distance perpendicular 
to the columns, such that the particles in these layers fit 
into the gaps in the layers below or above. In this way 
a higher density than Eq. ( |18[ ) is achieved. The columns 
of the third crystal phase (B) resemble braids with alter- 
nating tilt direction of the particles within each column. 
Because of this tilt Si and S2 have values between and 
1, that depend on D. Furthermore, the inverted braids 
structure (IB), that has < S2 < 1 and Si = 0, can be 
obtained by flipping one half of the columns of the braid- 
like phase (B) upside down. These braid-like columns 
piece together in such a way that the particles are inter- 
digitated. In other words, this phase is related to the B 
phase in exactly the same way as the IX phase is related 
to the X phase. Finally, in the paired face-centered-cubic 
(fcc^) phase, pairs of hemispheres form sphere-like ob- 
jects that can rotate freely and that are located at the 
lattice positions of an fee crystal. The density at close 
packing is 2^/2 /a^, i.e. twice the density of fee. 

In Fig. [8] the results of the pressure annealing method 
are shown, along with the known packing fraction of the 
perfect hexagonal columnar phase (col). Since the colum- 
nar phase has positional degrees of freedom and the fcc^ 
phase has rotational degrees of freedom, we expect these 
phases to have a higher entropy (lower free energy) than 
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FIG. 8. Packing diagram: maximum packing fraction (<j!>) of 
various crystal phases as a function of the thickness (D) of the 
bowls. The points are the results of the pressure annealing 
simulations. The thin dot-dashed lines are obtained from the 
pressure annealing results by slowly increasing or decreasing 
D as described in Sec. |IID[ except for the IX phase (thin 
dashed line with open squares) and the X phase (thin solid 
line with filled squares), for which the packing fraction can 
be expressed analytically. The thick lines denote the pack- 
ing fractions of the perfect hexagonal columnar phase (col) 
and the paired fee phase (fcc^). Any points that lie below 
these lines are expected to be thermodynamically unstable 
(see text). 



any crystal phase with the same or lower maximum pack- 
ing fraction whose degrees of freedom have all been frozen 
out. Therefore, any crystal structure with a packing frac- 
tion below the thick lines in Fig. [8] is most likely thermo- 
dynamically unstable. At first, we were unable to find the 
fcc^ using the pressure annealing method as described in 
Sec. |IID[ However, if we increase the pressure slowly to 
100fcsT/(T^ in simulations of 12 particles, we did observe 
the fcc^ phase for hemispherical particles {D = cr/2). In 
these simulations at finite pressure, it is important to 
constrain the length of all box vectors such that they 
remain larger than say 1.5a. Otherwise the box will be- 
come extremely elongated, such that the particles can 
interact primarily with their own images. When a parti- 
cle interacts with it is neighbors, the Gibbs free energy 
G = F + PV decreases, because the volume decreases 
without any decrease in entropy due to restricted trans- 
lational motion (if a particle moves, its image moves as 
well, so a particle translation will never cause overlap of 
the particle with its image). The decrease in Gibbs free 
energy is of course an extreme finite size effect, which 
should be avoided if we wish to predict the equilibrium 
phase behavior. For the pressure annealing simulations 
at very high pressures, these effects are not important, 
because the entropy term in the Gibbs free energy is small 
compared to PV. We did not attempt to find the colum- 
nar phase using the modified pressure annealing method, 
as we were only interested in finding candidate crystal 
structures. Furthermore, the columnar phase was already 
found in more standard simulations with a larger number 
of particles. 



C. Free energies 

In order to determine the regions of the stability of the 
fluid, the columnar phase and the six crystal phases, we 
calculated the free energies of all phases as explained in 
the Methods section. The results of the reference free 
energy calculations are shown in Tbls. |T]and|IT] 

We find that the columnar phase with all the parti- 
cles pointing in the same direction is more stable than 
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FIG. 9. Dimensionless free energies l3Fa^/V for hard bowls 
with L — 0.3(T and the fluid-columnar, columnar-IX and IX- 
IB coexistences, which were calculated using common tangent 
constructions. The columnar phases is denoted "col". The 
irrelevant free energy offset is defined in such a way that the 
free energy of the ideal gas reads I3F/V — p(log(po"^) — 1). 
The free energies of the various phases are so close, that they 
are almost indistinguishable. 



the inverted columnar phase, where half of the columns 
are upside down. However, the free energy difference be- 
tween the two phases is only 0.013±0.002A:bT per particle 
at = 0.5193 and D = O.Scr. Based on this small free 
energy difference we do not expect polar ordering to oc- 
cur spontaneously. Similar conclusions, based on direct 
simulations, were already drawn in Ref. [2]. 

The densely-packed crystal structures in Fig. [7] at 
D < 0.3, the worm-like fluid phase (Fig. [2]) and the 
columnar phase (Fig. ^ show striking similarity in the 
local structure: in all these phases the bowls are stacked 
on top of each other, such that (part of) one bowl fits 
into the dent of another bowl. As a result, the free ener- 
gies and pressures of the various phases, are often almost 
indistinguishable near coexistence. For this reason it was 
sometimes difficult to determine the coexistence densities 
for D < 0.3ct. Exemplary free energy curves for the var- 



phase 
fluid-col 



D/a 
1) 



/OfluidO" 

1.461 



4.679 



/diflF 

7.33272 



phases 


D/g 


0fluid 


0CO1 


/diff 


fluid-col 


0.1 


0.1780 


0.2848 


3.2630(7) 


fluid-col 


0.2 


0.3116 


0.4674 


3.268(2) 


fluid-col 


0.3 


0.3760 


0.5193 


3.802(1) 


fluid-inv col 


0.3 


0.3760 


0.5193 


3.8155(8) 


fluid-col 


0.4 


0.4440 


0.5772 


5.843 



TABLE I. Free energy differences, /diff = (-Fcoi(pcoi) — 
Ffinidip fluid) / {N ksT) , between the (inverted) columnar phase 
at density pcoi or packing fraction 4>coi and the fluid phase at 
Pffuid or (^fluid. In the column "phases", "col" denotes the 
columnar phase and inverted columnar phase is abbreviated 
to "inv col" . 



phase 


D/o 


4> 


/cxc 


IX 


0.3 


0.6669 


15.505(4) 


IB 


0.3 


0.6971 


18.407(3) 


IX 


0.4 


0.6177 


12.52(1) 


IB 


0.4 


0.6170 


13.195(2) 


IX 


0.45 


0.6768 


17.918(2) 


IB 


0.45 


0.6662 


14.9873(4) 


fcc^ 


0.45 


0.6192 


12.8591(5) 


IX' 


0.45 


0.6950 


18.170(5) 


fcc^ 


0.5 


0.5455 


8.7673(7) 


IX' 


0.5 


0.5597 


10.854(3) 



TABLE II. Excess free energies, /oxc = {F - Fa)/{NkBT), 
of the various crystal phases, where Fid is the ideal gas free 
energy. The various crystal phases are labeled as in Fig. [7] 



ious stable phases consisting of bowls with D 
shown in Fig. [9J 



3(7 are D /a phase 1 phase 2 p\(j^ p2<J'^ pPa^ 



D. Phase diagram 

In Fig. |10[ we show the phase diagram in the packing 
fraction </> - thickness D/a representation. The pack- 
ing fraction is defined as cj) = ^{a"^ - Da + ID^)N/V. 
For D/a < 0.3, we find an isotropic-to-columnar phase 
transition at intermediate densities, which resembles the 
phase diagram of thin hard discs ^^. However, the fluid- 
columnar-crystal triple point for discs is at a thickness- 
to-diameter ratio of about L/a ~ 0.2 — 0.3, while in our 
case the triple point is at about D/a ^ 0.3 — 0.4. The 
shape of the bowls stabilizes the columnar phase com- 
pared to the fluid and the crystal phase. We find four 
stable crystal phases IX, IB, IX' and fcc^, while we had 
six candidate crystals. The two phases that were not sta- 
ble are the X and B crystals, which are very similar to 
the stable IX and IB crystals respectively, except that 
X and B have considerable lower close packing densities. 
Therefore, one could have expected these phases to be 
unstable. On the other hand, we observe from the phase 
diagram, that IX is stable at intermediate densities for 
0.25ct < D < 0.45(7, while IB packs better than IX. In 
other words, stability can not be inferred from small dif- 
ferences in packing densities. 

Almost all coexistence densities were calculated by em- 
ploying the common tangent construction to the free en- 
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FIG. 10. Phase diagram in the packing fraction (</!>) versus 
thickness (D) representation. The light gray areas are coex- 
istence areas, while the state points in the dark gray area are 
inaccessible since they lie above the close packing line. IX, IB, 
IX' and fcc^ denote the crystals as shown in Fig. ItI "F" is the 
fluid and "col" is the columnar phase. The lines are a guide 
to the eye. Worm-like stacks were found in the area marked 
"worms" bounded from below by the dashed line. This line 
denotes the probability to find a particle in a cluster that 
consists of more than two particles, PstackC^i > 2) = 1/2. 
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TABLE III. Reduced densities, pressures and chemical poten- 
tials /i* = /3/i — ln(At Ar/cr'^) of the coexisting phases for hard 
bowl-shaped particles with thickness D. 



ergy curves, except for the col-IX coexistence at _D = 
O.Ict and 0.2cr. At these values of D the transition oc- 
curs at very high pressures, while the free energy of the 
columnar phase is calculated at the fluid-col transition, 
which occurs at a low pressure. To get a value for the free 
energy of the columnar phase we would have to integrate 
the equation of state up to these high pressures, accu- 
mulating integration errors. Furthermore, we expect the 
coexistence to be rather thin, which would further com- 
plicate the calculation. So, instead we just ran long vari- 
able box shape NPT simulations to see at which pressure 
the IX phase melts into the inverted columnar phase. As 
the free energy difference between the inverted colum- 
nar phase and the columnar phase is small, we assume 
that this is the coexistence pressure for the col-IX tran- 
sition, although technically it is only a lower bound. The 
density of the columnar phase at this pressure is deter- 
mined using a local fit of the equation of state. All co- 
existences are tabulated in Tbl. IIIII We draw a tentative 
line in the phase diagram to mark the transition from a 
structureless fluid to a worm- like fluid i.e. a fluid with 
many stacks. In a dense but structureless fluid, stacks 
of size 2 are quite probable, but larger stacks occur far 
less frequently. We calculate the probability to find a 
particle in a stack that contains more than 2 particles 
"Pstackl"- > 2) = 1 - T'stackll) - ^stack(2) and define the 
worm-like phase by the criterion 7^stack('^ > 2) > 1/2 
in Fig. [To] We do not imply that the transition to the 
worm-like phase is a true thermodynamic phase transi- 
tion; the transition is rather continuous. The type of 
stacks in the fiuid changes from worm-like for D — 0.3ct 
to something resembling the columns in the braid-like 
crystals B and IB (see Fig. [?]) for £» = 0.4a. There- 
fore, the region of stability worm-like phase was ended at 
D = 0.35(7, where there are similar amounts of braid-like 
and worm-like stacks. 
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IV. SUMMARY AND DISCUSSION 

We have studied the phase behavior of hard bowls in 
Monte Carlo simulations. We find that the bowls have a 
strong tendency to form stacks, but the stacks are bent 
and not aligned. We measured the equation of state and 
the compressibility in Monte Carlo NPT simulations. 
The pressure we obtained from these simulations is con- 
cave for some range of densities for deep bowls. This is 
due to the increase in free volume when large stacks form. 
Using gc{z), the pair correlation function along the direc- 
tion vector, we showed that the concavity of the pressure 
coincides with a dramatic change in structure from a ho- 
mogeneous fluid to the worm-like fluid. We measured 
the three-dimensional stack length distribution in the 
simulations. When the pressure is increased slowly, the 
deep bowls spontaneously order into a columnar phase 
in our simulations. This poses severe restrictions on the 
thickness of future bowl-like mesogens (molecular or col- 
loidal) , which are designed to easily order into a globally 
aligned lyotropic columnar phase. We determined the 
phase diagram using free energy calculations for a parti- 
cle shape ranging from an infinitely thin bowl to a solid 
hemisphere. We find that the columnar phase is stable 
for D < 0.3a at intermediate packing fractions. In ad- 
dition, we show using free energy calculations that the 
stable columnar phase possesses polar order. However, 
the free energy penalty for hipping columns upside down 
is very small, which makes it hard to achieve complete 
polar ordering in a spontaneously formed columnar phase 
of bowls. 
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of the bowl D, in the following way: 




f?n f?, 


(Al) 


"^^ "^^ ' 2iR,-D) 


02 = arcsin(i?i/i?2) 


(A2) 


L = i?2Cos(e'2). 


(A3) 



Overlap occurs if either of the two parts of the sur- 
face of a bowl overlaps with either of the two parts of 
another bowl. So we have to check four pairs of infinitely 
thin (and not necessarily hemispherical) bowls, labeled 
i and j, for overlap. The algorithm for two such sur- 
faces that are equal in shape was already implemented 
by He and Siders jSlj as part of their overlap algorithm 
for their "UFO" particles, which are defined as the in- 
tersection between two spheres. An equivalent overlap 
algorithm was used by Cinacchi and Duijneveldt [12^ to 
simulate infinitely thin contact lense-like particles, but 
the overlap algorithm was not described explicitly. We 
can not use one of these algorithms, since the two parts of 
the surface of our particle are unequal in shape. There- 
fore, we implemented a slightly different version of the 
overlap algorithm, which we describe in the remainder of 
this section. In our overlap algorithm, the existence of a 
overlap or intersection between two infinitely thin bowls 
is checked in three steps. 

• First, we check whether the full surfaces of the 
spheres intersect, i.e. \Ri — Rj\ < rij = |rj — r^j < 
Ri + Rj. If this intersection does not exist, there is 
no overlap, otherwise we proceed to the next step. 

• Secondly, we determine the intersection of the sur- 
face of each sphere with the other bowl. The inter- 
section of bowl i with the sphere of bowl j exists 
if 



kij +C(t)ij\ < 



for C = 1 or — 1, where 



m 



cos I 



^1 



zr^j -Ttj 



and 



cos{ujij) 



Ui 



(A4) 

(A5) 
(A6) 



Appendix A: Overlap algorithm 

The overlap algorithm for our bowls checks whether 
the surfaces of two bowls intersect. Fig. [T] shows that 
the surface of the bowl consists of two parts. Part p of 
the surface contains the part of the surface of the sphere 
of radius Rp, within an angle 9p from the z-axis, where 
p = 1 denotes the smaller sphere and the larger sphere is 
labeled p = 2. We set 9i = tt/2, to get a hemispherical 
outer surface. The edges of both surfaces have to coin- 
cide, such that our particles have a closed surface. Using 
this restriction L, 02 and i?2 can all be expressed in terms 
of the radius of the smaller sphere, i?i, and the thickness 



see Fig. [TT^ . This intersection is an arc, which is 
part of the circle that is the intersection between 
the two spheres. If in fact this arc is a full circle and 
the other particle has a nonzero intersection, the 
particles overlap. This is the case when Eq. (A4) 
holds for C = 1 o,nd, ^ = — 1. If, on the contrary, 
either of the two arcs does not exist, there is no 
overlap. Otherwise, if both arcs exist, but neither 
of them is a full circle, proceed to the next step. 

• Finally, if the two arcs overlap there is overlap, oth- 
erwise the particles do not overlap. The arcs over- 
lap if 



m 



<h\ + \ijl 



(A7) 
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where 



COs(q;,j) = 
cos(7i) = 



n 



_L 



n" 



-L||„-L 



cos 



ni 



^) 



COS! 



,-)cos(wy) 



sm{(j)ij) sm{ujij) 



(A8) 
(A9) 



where n^ = n^ — (ry-ni)rij/r?- and the expressions 
for 7j and nj^ are equal to the expressions for 7^ and 
n^ with i and j interchanged. The arcs together 
with the relevant angles are drawn in Fig. [lib. 



The inequalities (A4) and (A7) are expressed in cosines 



and sines using some simple trigonometry. In this way no 
inverse cosines need to be calculated during the overlap 
algorithm. 

For D — 0.5cr the bottom surface is a disk rather than 
an infinitely thin bowl. So the overlap check consists 
of bowl-bowl, bowl-disc and disc-disc overlap checks. 
For brevity, we will not write down the bowl-disk over- 
lap algorithm, but it can be implemented in a similar 
way as the algorithm for bowl-bowl overlap described 



above. The disk-disk overlap algorithm was already im- 
plemented by Eppenga and Frenkel |32j . 



(b) in 





FIG. 11. The relevant lengths and angles which are used in 
the first and second steps (a) and in the third step (b) of the 
overlap algorithm. Shown are bowl i and (part of) the sphere 
of bowl j (a), the arcs of i and j and the circular intersection 
of the spheres of i and j (b). In (a) r^j lies in the plane, while 
the plane of view in (b) is perpendicular to Tij. In this case, 
the sphere of particle j overlaps with bowl i, but the arcs do 
not overlap, so particle i and particle j do not overlap. 
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